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ABSTRACT 

A non-resonant instability for the ampli ficati o n of t he interstellar magnetic field in young Supernova 
Remnant (SNR) shocks was predicted bv iBelll (|2004f ). and is thought to be relevant for the accelera- 
tion of cosmic ray (CR) particles. For this instability, the CRs streaming ahead of SNR shock fronts 
drive electromagnetic waves with wavelengths much shorter than the typical CR Larmor radius, by 
inducing a current parallel to the background magnetic field. We explore the nonlinear regime of 
the non-resonant mode using Particle-in-Cell (PIC) hybrid simulations, with kinetic ions and fluid 
electrons, and analyze the saturation mechanism for realistic CR and background plasma parameters. 
In the linear regime, the observed growth rates and wavelengths match the theoretical predictions; 
the nonlinear stage of the instability shows a strong reaction of both the background plasma and the 
CR particles, with the saturation level of the magnetic field varying with the CR parameters. The 
simulations with CR-to-background density ratios of ncii/^b = 10 -5 reveal the highest magnetic field 
saturation levels, with energy also being transferred to the background plasma and to the perpendic- 
ular velocity components of the CR particles. The results show that amplification factors > 10 for 
the magnetic field can be achieved, and suggest that this instability is important for the generation 
of magnetic field turbulence, and for the acceleration of CR particles. 
Subject headings: cosmic rays — instabilities — magnetic fields — supernova remnants 



1. INTRODUCTION 

Very energetic CRs (~ 10 14 eV to 10 15 eV) are 
thought to be accelerated in SNR shocks. There is di- 
rect evidence that electrons are accelerated up to en- 
ergies of 10 14 eV at SNR sites i dKovama et all 119951: 
Allen et a l. 1997; Tanimo ri et al.lfl998l; lAharonian et al.l 
2QQl[lNaito et al. 1999; Aharonian 199*1 iBerezhko et al.l 
2003; IVink fc Lamingl [2003h . and the measured power 
law spectra of the CRs indicates Diffusive Shock Accel- 
eration (DSA) as the most likely mechanism responsi- 
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ble for the acceleration (lAxford et al.l 119771 : [Belli 119781 : 
iBlandford fc O striker! 1 1978[ ). The acceleration of these 
particles up to energies of ~ 10 15 eV through the DSA 
mechanism requires the existence of magnetic fields much 
stronger than the typical Bq ~ 3 fiG. These strong 
fields have also r ecently been inferred fr om observations 
(lLongairl Il994t IBerezhko et~aTI 12003b IVink fc Lamingl 
l2QQ3h : their existence, along with the requirement of 
stronger fields for the DSA mechanism, suggests that a 
magnetic field amplification mechanism is in operation. 

A possible amplification mechanis m through a non- 
resonant instability was suggested in IBelll (I2QQ4D. follow- 
i ng pr evious work on the resonant mode in lLucek fc Belli 
(|2000l ), a nd la t er ext ended to include multidimensional 
effects in Bell| (|2QQ5h . The non-resonant streaming in- 
stability, part of a class of stre aming instabilities de- 
rived in IWinske fc Lerovl (|l984f ). can be described by 
considering a MHD model for the background plasma, 
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and a CR induced current imposed externally; the feed- 
back of the electromagnetic fields on the CR particles 
is thus neglected. Although in the linear stage of the 
instability A max <C tlcr (A m ax the fastest growing wave- 
length and tlcr the typical Lar mor radius of th e CRs) , 
recent full PIC simula t ions by [Niemiec et al.l (12 008). 
Riauelme fc S pitkovskvl (|2009f ). lOhira et all (|2009l ) and 
St r oman et al.l ((2009) indicate that the feedback mecha- 
nism of the fields on the CRs is important in the non- 
linear stage, and suggest that a careful study of the in- 
stability in this regime is important to determine the 
saturation levels of the ma gnetic field. Rec ent works by 
lAmato fc BlasH ()2QQ9h and iLuo fc Melrosd (|2009D . using 
kinetic theory, have also shown how the saturation of this 
non-resonant mode depends on the details of the particle 
distributions. 

The analysis of the behavior of the non-linear stage 
of the instability is very complex; full assessment of the 
saturation level of the magnetic field would imply a first- 
principle calculation of the shock formation, and the long 
term evolution of the shock precursor, including the self 
consistent acceleration of particles. The saturatio n of 
the in stability was thus first assessed numerically in lBelll 
(2004), where an external current was used to model the 
CRs, and an MHD model was used to simulate the back- 
ground plasma. The saturation mechanism found was 
due to the tension of the magnetic field, which grows 
faster than the driving term; also, the size of the sim- 
ulation b ox was see n to limit further growth of the in- 
stability (Bell 20 04]). More recent PIC simulation re- 
sults in Niemiec et al.l ((2008) actually show a satura- 
tion level of SB/B ~ 1 5 but consider parameters such 
that 7max/^ci <^ 1 (the ratio of the growth rate of the 
fastest growing mode to the CR cyclotron frequency) is 
not strictly maintained, which is a requirement for the 
development of the non-r esonant parallel mode. The ful l 
PIC simulation results in IRiquelme fc Spitkovskvl (|2009l ) 
and IStroman et al.l (j2QQ9h , taking into account the feed- 
back of the electromagnetic fields on the CR particles, 
and using rrii/m e ratios up to 100 and ncn/^h down 
to 4 x 10 -3 , imply a saturation level of SB/B ~ 10, 
which occurs when the relative drift between the CRs and 
the b ackground plasma decreases. Also, in lOhira et"aTI 
(2009), similar saturation levels for the magnetic field are 
found. 

For PIC simulations, therefore, the instability satu- 
rates when CRs loose some of their bulk momentum 
to the background plasma, which starts to drift more 
rapidly in the direction of CRs until the velocities of 
the two species converge and the driver for the insta- 
bility is elimin a ted. As shown in kinetic modeling by 
ILuo fc Melrosej (|2009[ ). the reduction of the CR stream- 
ing motion is due to CR resonant diffusion in non- 
resonant mag netic turbulence, wh ich has been recently 
observed in Stroman et all (|2009l ). The nonlinear evo- 
lution of the turbulence observed in fully kinetic simu- 
lations is also in qualitative agreement with the predic- 
tions of quasi-linear calculations for no n-resonant modes 
derive d for non-relativistic beams in IWinske fc Lerovl 
(|l984f ). However, the accurate modeling of the late time 
evolution of the system, beyond the saturation, is lim- 
ited in the kinetic simulations, when the assumption of 
plasma homogene ity is no longer valid. As argued in 
(|Ohira et al.ll2009l ). when the background plasma veloc- 



ity approaches the CR bulk velocity, the density in a real 
scenario also changes, modifying the underlying condi- 
tions for the simulation model. This reinforces the ne- 
cessity for a thorough analysis of the behavior of this 
instability in its non-linear stage, covering a wide range 
of dynamical scales, and exploring the saturation limits 
in different regimes. 

Here, we present multi-dimensional simulation results 
of the non-resonant streaming instability, using the ki- 
netic ion flu id electron hybrid m odel implemented in 
dHybrid (see iGargate et al.l ((2007) for numerical imple- 
mentation details). An important advantage of the hy- 
brid simulations is to enable the study of the instabil- 
ity on the ion time scale, neglecting the high-frequency 
modes associated with the electrons; realistic density ra- 
tios down to ncR/nb = 10 -5 can then be used, along 
with realistic ion-to-electron mass ratios, and simula- 
tions can be run to a point well beyond the linear stage, 
into the saturated state. The variation of the den- 
sity ratio has an impact on the behavior of the non- 
resonant mode, with lar ger v alues implying the gener- 
ation of different modes (iBehl 120041: iNiemiec et alJ l2QQ8t 
IRiquelme fc Spitkovskvl l2QQ9h . Though the saturation 
mechanism is indepen dent of the initial linear instability 
([Stroman et al.l l2009). these modes lead to smaller satu- 
ration amplitudes and thus may not be directly relevant 
for the typical SNR shock scenarios. Also, in the hybrid 
simulation results shown, both the background and the 
CR ions are modeled kinetically, which enables the study 
of the feedback mechanism on the CR population. Re- 
sults can then be easily compared with both the external 
current-driven MHD simulations, where feedback on the 
CR population is not modeled, and the fully kinetic sim- 
ulations, which use small ion-to-electron mass ratios, and 
larger ucr/ti^ ratios. 

This paper is organized as follows. In section 2, the pa- 
rameters used in the simulations are discussed, and the 
results concerning the evolution of the linear and nonlin- 
ear stages of the instability are presented. The satura- 
tion mechanism is discussed in section 3, and compared 
with the most recent results in the literature. Finally, in 
the last section, we present the conclusions and outline 
future research directions. 

2. EVOLUTION OF THE NON-RESONANT INSTABILITY 

The original theoretical model of the non-re sonant 
streaming instability, as detailed in Bell (|2004f ). con- 
siders that a magnetized background plasma with den- 
sity nb is stationary in the SNR shock upstream refer- 
ence frame. The CR particles stream along the back- 
ground magnetic field lines B\\ with a drift velocity v s h, 
similar to the shock velocity, generating a current that 
drives the growth of the perpendicular components of 
the magnetic field B \ . For the canonical parameters 
in the literature (jBelll 120041 ). of B\\ = 3/iG, n b = 
lcm -3 , v s h = 10000 km/s, and tlcr ~ 1-1 x 10 13 m, 
the Alfven velocity is = 6.6 km/s, and the Alfven 
Mach number is Ma = v s h/vA = 1515 ^> 1. Un- 
der these conditions, the reduced wave equation uo 2 — 
v\k 2 ± |£?||j||| (n^m) = 0, where m is the mass of 
the background ions, j\\ is the zeroth order current im- 
posed by the CRs, k is the wavenumber and uo is the fre- 
quency of the mode, yields a purely growing mode with a 
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Fig. 1. — Perpendicular magnetic field component B z (frame a 
and b), and the background plasma density (frame c and d) for 
run &2- Frames a) and c) correspond to time 7.54 r y~ 1 in the linear 
stage; frames b) and d) to time 16.21 7 — 1 in the nonlinear stage. 



TABLE 1 

Key simulation parameters 



Run 


Dim. 






^iso 


A 


2D 








Bi 


2D 


io- 3 


157 v A 


800 v A 


B 2 


2D 


io- 4 


1570 v A 


800 v A 


B 3 


2D 


icr 5 


15700 v A 


800 v A 


Ci 


2D 


icr 3 


157 v A 


160 v A 


c 2 


2D 


io -3 


157 v A 


320 v A 


Cs 


2D 


icr 3 


157 v A 


500 v A 


V 


3D 


io- 3 


157 v A 


800 v A 



1/2 

growth rate 7 = [Bujuk/ (rib m) — v a^ 2 ] •> an d a maxi- 
mum growing wavenumber fc max = 1/2 \B\\ j\\ | / (rib m v\) 
with a growthrate of 7 max = ^max^A- For the non- 
resonant mode to be predominant, then 1 < Atlcr < 
C v sh/ V h witn C = \B\\j\\\r LC R/(nbmv^ h ); for the canon- 
ical parameters it follows that /c m ax ^lcr = 3691, and 
C v sh/ V A ~ 7381. Moreover, 7max/^ci <C 1 is maintained, 
which is an essential requirement for the development of 
the non-resonant parallel mode. 

The results are presented in normalized units, with 
the magnetic field normalized to the background field 
Bq = Bn, the density normalized to the background 
plasma density no = nt>, the velocities normalized to 
the Alfven velocity va — Bo/(nom p liq) 1 ^ 2 ? the spatial 
dimensions normalized to the ion inertial length c/uj p i, 
and time normalized to the inverse maximum growth 
rate 7m ax - The normalized background magnetic field in 
the simulations is Bn = B x = lBo, and the normalized 
background plasma density is rib — l n o, corresponding 
to SNR scenarios where B\\/rib = 0.4/iGcm 3 . 

The main parameters for the simulations are described 
in Table [U The CRs are modeled through a constant 
external current in run A, and with kinetic particles 
in all other runs (the kinetic-CR runs). The value 
of the initial CR current is always j = uqr e v s h = 
7r/20 [no e va] (where [no e va] are the normalized current 
density units), which yields a maximum growing wave- 
length A max = 80c/cj pi (/c max = 7r/40cj pi /c). For the 
kinetic-CR runs, the total velocity in the CR population 



is initially set to v s h + Vi SO , with V[ so the isotropic velocity 
component, so that the whole CR population drifts with 
the shock velocity, as observed in the upstream reference 
frame. For runs B\ through 63 the CR-to-backghround 
density ratio was varied, along with the drift velocity v s h, 
and for runs C\ through C3, the V{ so was varied, enabling 
the analysis of the dependence of the instability on the 
scale separation between the typical Larmor radius of the 
CRs, and A max . Finally, in run T>, the parameters from 
run B\ were used in a 3D setup. For all the runs the time 
step was chosen to be At = 0.001 co^ 1 , as to account for 
the cavitation in the background density plasma, since 
low density zone s can render the hybrid (and MHD) algo- 
rithms unstable (|Lipatovll2002t iGargate et al.ll2007l i. Full 
convergence tests have been performed to assure that the 
saturation level of the magnetic field was not affected by 
the choice of the time step. 

The simulation box size is L x = L y = 1280c/cc;pi for 
the 2D runs, and L z = 640c/cj p i for the third dimension 
in the 3D run presented. The L x and L y dimensions are 
then ~ 16 x A max , corresponding to L x /L^ ~ 1, where 
Ldiff ~ D/v^ is the CR diffusion length scale in the 
Bohm diffusion limit for the typical CR velocities used 
in the simulations. This choice of parameters then al- 
lows for the CR Larmor radius to be fully resolved, and 
is such that the non-resonant mode wavelength can grow 
to be of the order of tlcr in the non-linear regime, en- 
abling the analysis of the back-reaction of the generated 
electromagnetic fields on the CR particles. 

The diverse parameters chosen for the four sets of runs 
(A through V) does not directly reflect typical physical 
conditions found in young SNR shock precursors for all 
the runs; instead a wide dynamical range of parameters 
was analyzed. The results of ru n A are mo st directly 
comparable to the simulations by Bel| (|2004f h in which 
the CR current was imposed externally, thus excluding 
any feedback of magnetic turbulence on theCRs. For the 
B set of runs, the driving current is maintained constant 
while the driving energy of the CR beam is increased 
from run B\ to run B3 by increasing v s h and decreasing 
ncR; as such ^i SO /^sh is smaller than the typical values 
found in young SNR shock precursors, especially for runs 
B2 and S3. The C set of runs explores the behavior of 
the instability when progressing from the non-resonant 
to the resonant limit of the CR streaming instability by 
lowering the Vi SO velocity and thus the nLCR/Amax ratio. 

Comparison of results between 3D and 2D runs showed 
no significant difference in the development of the insta- 
bility. This motivated the analysis of the system using 
a 2D simulation setup. For run V, the measured wave- 
length of the instability was A = 80 c/cj p i, and the growth 
rate 7 ~ 0.889 of the non-resonant mode is the same as 
for all the runs in the B set. Figure [T] shows the evolu- 
tion of the B z magnetic field component and the density 
of the background plasma for run B2] Fig. Q] a) and 
c) correspond to the linear stage and show the well de- 
fined wavelength A = 80c/o;pi in the magnetic field, and 
a low amplitude modulation of the background density. 
In the nonlinear stage, cavities are formed in the back- 
ground plasma, due to turbulent plasma motions induced 
by the electromagnetic waves being generated. The well 
defined wave structure of the linear stage disappears in 
Fig. Q]b), with magnetic field compression zones being 
formed, which are correlated with the density enhance- 
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Fig. 2. — Evolution of the perpendicular magnetic field com- 
ponent By in time for run &2- Frame a) shows lineouts of the 
power distribution in k space (frame b) at times: 27 -1 (black line), 
67 -1 (blue line), 87 -1 (green line), and 16 7" 1 (red line). Frame 
d) shows the growth rates for k x = 0.078 w p i/c (solid line), k x = 
0.044 u; pi /c (dashed line), and for the 0.029 < k x < 0.132 o; pi /c 
high-growth band (dotted line). Frame b) is obtained by fourier 
transforming frame c) in the spatial dimension. 

merits at t he same positions, Fig. Ud), confirming the 
results of ([Riquelme fc Spitkovsky|[2QQ9f ) for lower mass 
ratios. Also, very similar results are observed for runs B\ 
and S3. 

In Fig. [2j the main characteristics of the instability 
that agree with expectations of a quasi-linear theory are 
observable. At t ~ 97 -1 the instability becomes non- 
linear, and saturation of the linear mode is observed a 
short time after, around t = 10 7 -1 (Fig. [2]b) and d). 
In the linear stage, the growth rate from Fig. [2]d) yields 
7 = .889 which deviates from the theoretical value by 
11%. It is also clear from Fig. [2] a) that the peak power 
is for k ~ 0.078 cjpi/c up to t = 87 -1 (dotted line), and 
that at t = 16 7 -1 the peak shifts to fc ~ 0.044 cj p i/c 
(dashed line). 

In the nonlinear stage, Fig. [2]b) and c), the increase in 
power over large fc is consistent with the magnetic field 
structures in Fig. [T]b), and indicates an increase in the 
turbulence level in B z . Figure [2] c) also shows that the 
waves move in the x direction at a velocity v ~ 1.8 va in 
the nonlinear stage. 

3. SATURATION MECHANISM 

The maximum energy level in the perpendicular mag- 
netic field components, at the end of the linear stage, is 
similar for the runs Bi, B2 (see Fig. [3]), and B3. In the 
linear stage of the instability, the behavior of run A is 
identical to the kinetically-driven runs, B\ and B2. In 
the nonlinear stage, however, the magnetic field energy 
peaks at t = 10 7" 1 for runs B\ and B2, while for run 
A the magnetic field energy still grows for t > 10 7" 1 . 
The growth beyond t = 10 7" 1 for run A is due to the 
continuous injection of energy into the system; a similar 
growth is observed for run B3, at a slower rate, which 
is due to the excess of free energy in the CR popula- 
tion for this run. Finally, run C\ (Fig. [3] b) represents 
an hybrid scenario, between the resonant and the non- 
resonant modes of the instability (nLRcMmax ~ 1.63), 
and thus the magnetic field growth rate is smaller in the 
linear stage. 



10 u 

^ 10- 2 

LU 

1 10- 4 

sz 
LU 

10- 6 

10° 

T= ID" 2 
LU 

I 10" 4 

c 
LU 

10- 6 





b 























5 10 
Time [y" 1 ] 



5 10 
Time [y 1 ] 



Fig. 3. — Time evolution of the energy for the magnetic field 
(red lines), the CRs (black lines), and the background plasma (blue 
lines); dashed lines represent the component parallel to the back- 
ground magnetic field x, and the solid lines represent the perpen- 
dicular component. Frame a) shows the evolution for run A, frame 
b) shows the evolution for run C\ , frame c) shows the evolution for 
run Bi, and frame d) shows the evolution for run £>2- The energy- 
scale is normalized to the total energy in run B3. 
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Fig. 4. — Time evolution of the energy in the perpendicular com- 
ponents of the magnetic field (top frames), and the time evolution 
of the current density for the ion component of the CRs, and the 
ion component of the background plasma (lower frames). Frame 

c) shows the current density components for run Bi, and frame 

d) shows the current density components for run B2] the vertical 
arrows, from left to right, indicate points where the background 
plasma current, J x is: at a minimum, zero, at a maximum, and 
zero, respectively. Lines are labeled in the plot for ease of interpre- 
tation. The energy scale is normalized to the total energy in run 
B 3 - 

The CR population velocity distribution is nearly 
isotropic at the end of all the simulations (as seen in 
the upstream / simulation frame). In the linear stage, 
the growth rates are identical between runs in the non- 
resonant regime, and the background plasma also gains 
energy at the same rate in the perpendicular direction 
(Fig. [3]c) and d). Unlike what is predicted from theory, 
v x of the background plasma increases, and the parallel 
energy becomes comparable to the perpendicular energy. 
With this energy increase, the plasma is no longer mag- 
netized, and thus the nature of the non-resonant mode 
changes. 

In order to understand the dependence of the non- 
resonant mode on the scale-separation between A max 
and tlcr, the isotropic velocity was varied in runs C\ 
through C3. The typical CR Larmor radius varies as 
r LC R = {130.6, 261.3, 408.2} c/cjpi (and 653.2 c/u; pi for 
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run Bi). Figure [3] and Fig. [4] a) and b) show that for 
^LCR/Amax ratios greater than ~ 3.2 (run C2), the behav- 
ior is consistent with the non-resonant mode, and only 
run Ci (rLCR/Amax = 1.63) results in an hybrid scenario. 
Analysis of Fig. |4] c) and d) also shows that for runs B\ 
and 62, at t ~ 9 7 _1 , the J x current density component 
for the CRs is not close to zero. This means that the 
saturation mechanism does not depend directly on the 
amount of free energy in the CR population. Instead, 
the saturation of the instability is due to the energy gain 
of the background plasma, which results in an average 
bulk flow in the CR propagation direction, and in the 
de-magnetization of the plasma. 

The relation between J x for the background plasma, 
and both the onset and the peak in the magnetic field en- 
ergy level of the instability is highlighted in Fig. |4]c) and 
d). The peak of J x , for the background plasma, marks 
the end of the linear stage of the instability. Likewise, 
when J x is at a minimum (t ~ 37 -1 ), the instability 
is entering the initial linear stage. Finally, for run Bi, 
Fig. 2]a) and c), the energy in the CRs at the saturation 
point is transferred to the background plasma and to the 
magnetic field; for run B2, Fig. |4]b) and d), some of the 
energy is converted into perpendicular energy in the CR 
population, as J±_ increases from t ~ 10.8 7 -1 onwards. 

4. DISCUSSION AND CONCLUSIONS 

The identification of the nonlinear saturation mech- 
anism for the non-resonant streaming instability is im- 
portant for the determination of the maximum magnetic 
field amplification in SNR shock scenarios. Amplifica- 
tion factors of B±/B\\ ~ 10, with local temporary peaks 
of B±/B\\ ~ 25 are observable in the current hybrid 
simulations; similar result s were also recently shown in 
the kinetic s i mulat ions of Riquelme & Spitkovskyl (|2009f ) 
lOhira et al.l (|2009h and iStroman et al.l ([2009). with re- 
duced temporal and spatial scales, mass ratios and den- 
sity ratios. The hybrid simulation results presented ex- 
pand the current knowledge of the non-resonant stream- 
ing instability by extending the dynamical range under 
study, and by performing a detailed study of the insta- 
bility under different driving regimes. 

Leveraging on the hybrid simulations presented here, 
it is possible to scan the parameter space and analyze the 
nonlinear stage of the instability in detail, using realistic 
^CR/^b and m/m e ratios. Close examination of Fig. [3] 
and Fig. 2] shows a number of important points, relating 
to the saturation mechanism. As the ncn/n^ ratio is de- 
creased towards 10 -5 , the energy in the CRs increases, 
as in our setup the current is initialized at the same level 
in all simulations; the increase in the free energy of the 
CRs, associated with the higher drift velocity v s h, does 
not change the peak energy level of the magnetic field. 
The excess free energy is instead transferred to the back- 
ground plasma, and also to the perpendicular velocity 
components of the CR population (Fig. [3]b), run C\). 

The reaction of the background plasma is critical in 
the nonlinear stage of the instability. The average local 
Alfven velocity (va Bn/y/n calculated in each cell and 
averaged over the entire simulation box) is approximately 
constant over the linear development of the instability. 
In the nonlinear stage, va increases with the local mag- 
netic field B\\, and with the background density varia- 



tions (which are in phase), and thus va oc ^5B/B\\. The 
background plasma accelerates in the CR propagation di- 
rection, and the current density peaks, marking the end 
of the linear stage (Fig. 2J. At this point in time, the CRs 
v x velocity component is still v x ~ 1500 va, well above 
the Alfven velocity and above the background plasma 
drift velocity, which is sub-Alfvenic. This shows that the 
saturation is not dependent on the amount of free energy 
in the CRs, but instead results from the reduction of the 
scale separation of the background plasma, and the CRs. 

The CR's J± increases for run B2, Fig. |4]d), as it is 
substantially lower than J x initially, and the distribution 
becomes isotropic after the saturation, at t ~ 15 7 -1 . For 
run Bi in Fig. |4]c), the distribution is also isotropic at the 
end, but the CRs J± velocity component is not affected. 
Increasing the free energy in the CRs does not affect the 
saturation level of the magnetic field significantly, as long 
as the driving current is maintained. 

The hybrid simulation results presented thus show 
that, for a broad range of parameters, the nonlinear 
growth of the streaming instability is independent on 
the energy of the driving CR population, and depends 
only on the current carried by the CRs. The ampli- 
fication factor of ~ 10 for the perpendicular magnetic 
field above the seed Bn feed indicates that the mech- 
anism is relevant for the magnetic field amplification 
by CR particles in SNR shocks. Beyond the satura- 
tion point, for t > 16 7 -1 , the CR particle distribution 
becomes isotropic, with the instability saturating after 
10 e- foldings. The de-magnetization of the background 
plasma, as it gains energy, hinders further growth of the 
non-resonant mode; when both the CRs and the back- 
ground plasma are unmagnetized, the i nstability shou ld 
behave more like the Weibel instability (Weibel 1959). 

Our results thus indicate that the instability begins 
to saturate when the current carried by the background 
ions is similar to the current carried by the CR popula- 
tion (see Fig. [4j). This results in a final average velocity 
for the background plasma of ~ ncR/^b^sh so that 
Vb ~ 10 _5 v s h for realistic parameters. At this point, the 
instability is in the non-linear stage, and the magnetic 
field growth rate is much lower (t > 97 -1 in Fig. [3] c) 
and Fig [4] a) and c) for run Bi). In fact, even in the linear 
stage of the instability, energy is being transferred to the 
perpendicular components of the background plasma at 
a greater rate than into the magnetic field perpendicu- 
lar components (see Fig. |3]c), and thus at some point 
the background plasma de-magnetizes and the instability 
saturates. 

Previous simula ti on PIC r esults 

(iRiauelme fc Spitkoyskvl l2QQ9t lOhira et al.l 120091 : 
IStroman et all 12009) show similar saturation levels for 
the magnetic field 5B/B ~ 10: the instability saturates 
when the background plasma is moving with a bulk 
velocity comparable to the CR population drift speed, 
which is equivalent to saying that the currents carried by 
the background plasma and the CRs are similar (since 
ncR ~ n b , and th us Jcr = ncR, e ^ sh ~ Jb = n h ev h ). It 
has been claimed ( Oh ira et al.ll2009l ) that this saturation 
limit might be numerical rather than physical because 
in a real shock precursor a significant variation in the 
background plasma velocity should be concurrent with 
a variation in the background plasma density. However, 
here we show for ra$ ^> m e and ucr <C that the 
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instability saturates due to the de-magnetization of 
the background plasma; at that point the currents 
carried by the background ions and the CRs are 
comparable, although the velocities differ significantly 
(since in our hybrid simulations ncR <C implies that 
00 ^sh^CR/^b from the equality of the currents at 
the saturation). Our simulations thus reinforce the 
result SB/B ~ 10 at saturation, and indicate that the 
instability should be important in a young SNR shock 
precursor. The simulations present in the literature 
make assumptions about the physical scenarios, and 
therefore significant care should be taken when inter- 
preting and comparing the results. For instance, the 
back-reaction of the fields in the CRs is not accounted 
for when using external currents as drivers, in full PIC 
simulations mass ratios and density ratios far from 
the physical conditions have been explored, and the 
dynamics of electrons is not accounted for in hybrid 
simulations. Thus, our results further indicate that 
first-principle simulations of SNR shocks and particle 
acceleration should be attempted in a future work, 
since the different simulation methods can capture the 
instability, but differences might still be found when the 
simulation setup and the models are refined in order to 
more closely match the experimental conditions. 

Further analysis for SNR shocks will be possible by 
considering the saturation mechanism when unperturbed 
CRs are continuously injected into the simulation box. 



This will allow for an improved comparison of the SNR 
shock scenario, where CR particles are also thought to be 
continuously injected from the shock into the upstream 
medium. This will be explored in a future publication, 
leveraging on the unique characteristics of the hybrid 
model. Finally, with hybrid simulations, it will also be 
possible to determine the dependence of usual character- 
istics of the instability with the distance to the shock, 
and to study in detail the energy profile of the acceler- 
ated CR particles. 
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